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Abstract 

Carbon Capture and Storage (CCS) is a recently dis- 
cussed new technology, aimed at allowing an ongoing 
use of fossil fuels while preventing the produced CO2 
to be released to the atmosphere. CSS can be modeled 
with two components (water and CO 2) in two phases 
(liquid and CC>2)- To simulate the process, a multi- 
phase flow equation with equilibrium phase exchange is 
used. One of the big problems arising in two-phase two- 
component flow simulations is the disappearance of the 
nonwetting phase, which leads to a degeneration of the 
equations satisfied by the saturation. A standard choice 
of primary variables, which is the pressure of one phase 
and the saturation of the other phase, cannot be applied 
here. 

We developed a new approach using the pressure of 
the nonwetting phase and the capillary pressure as pri- 
mary variables. One important advantage of this ap- 
proach is the fact that we have only one set of primary 
variables that can be used for the biphasic as well as the 
monophasic case. We implemented this new choice of 
primary variables in the D UNE simulation framework 
and present numerical results for some test cases. 

1 Introduction 

In this work we address the mathematical modeling 
and numerical simulation of multiphase multicompo- 
nent flow in porous media with a special regard to CO2 
storage in geologic formations. Some people consider 
CO2 storage, e.g. in deep saline aquifers, as an im- 
portant factor in the effort to reduce the emission of 
greenhouse gases. Reliable simulation data is crucial 
for all stages of CCS projects. 



After the CO2 injection, several different trap- 
ping mechanisms lead to an entrapment of the CO2. 
Shortly after the injection, structural trapping through 
caprocks is the most important factor. Later solubil- 
ity trapping, where C0 2 is dissolved into water, and 
residual trapping get more important. After several 
thousand years, there could also occur mineral trap- 
ping caused by geochcmical reactions, but these are 
not considered in this work. 

The mathematical model describing CO2 injection in 
geologic reservoirs is a two-phase two-component flow 
in porous media, a system of coupled, nonlinear partial 
differential equations. We do not only have two differ- 
ent phases (liquid and CO2), but also two components 
(water and CO2) in each phase, as the solubility of 
the components in the phases has to be taken into ac- 
count. For an isothermal system we have to choose two 
primary variables and additional algebraic relations to 
close the system. 

A standard choice for the primary variables is the 
pressure of one phase and the saturation of the other 
phase. A great challenge in this context is the dis- 
appearance of the nonwetting phase, which has been 
studied in many recent papers, as the saturations can- 
not be used as primary variable here. A valid choice in 
the one-phase region would be one phase pressure and 
the solubility of CO2 in the liquid phase. 

Several Approaches to treat this problem exist, Class 
et al. (2002) switch primary variables depending on 
present phases, Jaffre et al. (2010) use complementarity 
conditions and Abadpour et al. (2009) extend the sat- 
uration to negative values. Bourgeat et al. (2010) use 
liquid phase pressure and water mass concentration as 
primary variables. 

In this study we present a new choice of primary 
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variables that is valid for the monophasic as well as 
the biphasic case and can easily handle the disappear- 
ance of one phase. One advantage of our approach is, 
that the chosen variables are continuous over material 
heterogeneities, if both phases are present. 

To simulate CSS, constitutive relations between 
physical properties like pressure and density are neces- 
sary. We summarize our choice of existing approaches. 

Numerical simulations for different test cases pre- 
sented in this work will show that this new approach 
handles various applications very well. We use our new 
approach to simulate CO2 injection into the subsurface. 
A recent benchmark from the MoMas group concen- 
trates on test cases arising from underground radioac- 
tive waste repository simulations. With our new set 
of primary variables we can also solve these problems 
efficiently 

2 Mathematical model of a isothermal two-phase 
two-component flow 

In this section we will consider a porous medium and 
derive a system of partial differential equations describ- 
ing two-phase two- component flow. For the sake of 
simplicity we use a constant temperature in this ar- 
ticle, but thermodynamic effects can be included into 
the model in a straightforward manner. We also as- 
sume that the salinity of the water is constant. 

2.1 Notation 

We have two phases a G {w, n}, wetting and nonwet- 
ting, and two components k € {a, b}, water and non- 
wetting component. 



p w , p n wetting and nonw. phase pressures 

S w , S n wetting and nonw. phase saturations 

Pmass.w, Pmass.n wetting and nonw. phase mass dens. 

Pmoi.w, Pmoi.n wetting and nonw. phase molar dens. 

/Lt w , ^ n wetting and nonw. phase viscosities 

x%, molar fraction of comp. in wet. phase 

x^, molar fraction of comp. in nonw. phase 

M a , M b molar mass of wet. and nonw. comp. 



where K is the absolute permeability, fc rw and k rn de- 
note the relative permeability functions and <? is the 
gravity vector. 

The phase saturations and molar fractions satisfy 

S n -\- 5\v 1, x w -\- x w 1, x n -\- x n 1. (3) 

The relation between the phase pressures is given 
through the capillary pressure by the Brooks-Corey or 
van Genuchten-Mualem model 

Pc{S w ) = p n - Pw (4) 

2.3 Diffusive flux 

Following Fick's Law, the diffusive flux of a component 
k in the phase a is given by 

3 a = _ -^pm,tt/ ) mol,oVj; a , (5) 

where D* m a is the diffusion coefficient of component 
k in phase a in a porous medium. 
Like pQ and [5] we assume 

3l+3% = (6) 

holds for simplicity, so we only need two diffusion co- 
efficients instead of four. 

2.4 Mass conservation 

Local equilibrium phase exchange of the components 
in the phases is assumed. Taking into account the con- 
servation of the amount of substance of each compo- 
nent and using Q, ([2| and ^ we get the following 
partial differential equations describing an isothermal 
two-phase two-component flow: 

-(- V • {pmol,w X w U w ~t- Pniol,nX n lin} 

+ V'{3l+3n} -<f = 0, 

,w ^w^w H~ Pmol,n ^n^n} 

+ V.{£+£} -«7 b -0, (7) 

where q a and q° are the source/sink terms for the com- 
ponents. 



2.2 Darcy's law 

The phase velocities u a are given by an extended 
Darcy's Law: 

U w = -K klW ^ W (Vj? w - jOmass.w " 9), (1) 
U n = -K k "^ Sl ^ (Vp n - Pmass.n ' ff), (2) 



3 Constitutive Relations 

We will now look at the special case of CCS where a 
liquid water phase and a liquid, gaseous or supercriti- 
cal CO2 phase are present. The components are water 
and CO2. In the following the different functions and 
Equations of State (EOS) to determine secondary pa- 
rameters are described. Additionally their dependence 
on other variables is given. 
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3.1 Solubility of components 



3.2 Densities 



The solubility of the components is influenced by the 
pressure p n and the temperature T of the system, the 
salinity s sa i of water also plays an important role: 



^* w (,Pn •, r ^; •''sal ) : ^ n (Pn ; T, 5 sa ] ) . 



(8) 



There exist different EOS for this system. We use the 
EOS by Spycher & Pruess [3], because in contrast to 
other models (for example, the EOS of Duan & Sun 
[I]) also the solubility of water in CO2 is described 
very well. Figure [T] and [2] show the solubility curves for 
different temperatures. 

The solubility of CO2 in the water phase increases 
fast with rising pressure up to the saturation pressure, 
above that it rises with a smaller rate. For tempera- 
tures below the critical temperature T cr it = 304.15 K, 
the state of the carbon dioxide changes from gaseous 
(below saturation pressure) to liquid which results in a 
not continuously differentiable sharp break at the tran- 
sition point. 
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Figure 1: Solubility x^, for different temperatures 
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Figure 2: Solubility i„ for different temperatures (s sa i = 0) 



For the density of the water phase the approach of Gar- 
cia [5] is applied. The density increases slightly for a 
larger fraction of CO2 in the water phase. The EOS 
of Duan [6] is used to calculate the density of the CO2 
phase, which strongly depends on the CO2 phase pres- 
sure, 

r (^¥i-^)i Pmass,n 



P. 



mass.w V 



Figure [3] shows the density of CO2 for different temper- 
atures. To convert mass density to molar density the 
phase composition has to be taken into account, 
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Figure 3: CO2 density for different temperatures 



3.3 Viscosities 

The viscosity of the water phase is computed with a 
function from Atkins [7], for the CO2 phase we use the 
approach of Fenghour & Vesovic |S]. Again the CO2 
phase viscosity strongly depends on the CO2 phase 
pressure, 

Figure [4] shows the viscosity of CO2 for different tem- 
peratures. 

3.4 Diffusion 

Following [5] , we use an approach suggested by Milling- 
ton & Quirk 



pm,Q 



(^) 10/ y 

/A2 «' 



for the diffusion coefficient in the porous medium, 
where describes the binary diffusion coefficient of 
component K in phase a. 
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Figure 4: CO2 viscosity for different temperatures 



4 Choice of primary variables 

A standard choice for the primary variables are one 
phase pressure and the saturation. In the one phase 
region (S n — 0), the system ^ degenerates to 

(f>df {Pmol,w £wl + V ' {Pmol.w ^w u w + j w } ~ (7* = °> 
4>9 t {pmol.w X w } + V ■ {pmolw 4 u w + jw} ~ <t = °- 

Using ([3]) and ([6l the system can be rewritten as a 
coupled groundwater-flow and transport problem 

<f>9t {Pmol.w} + V ' {/W,w«w} = <? a + <? b , (9) 

(f>d t {Pmol.w ^w} + V ' {Pmol,w ^ w u w + j w } ^ '■ 

With the disappearance of the nonwetting phase the 
saturation can no longer be used as primary variable 
and the standard choice of variables cannot be applied 
here. One natural set of variables for the one phase re- 
gion would be p w and x w . There are several approaches 



Primary Var. 



Method 



Sni Pn 
5w 5 Pw 5 



Extending the saturation 
to negative values (see jlOj). 
Using complementarity constraints 

(see[II]). 

p n , (S w or Xjj) Switching primary variables depen- 
ding on present phases (see [IS], [13]). 

Table 1: Several methods to deal with a disappearing 
nonwetting phase 



variables are used in the one phase and two phase re- 
gion, the variables are switched if a phase appears or 
disappears. 

Abadpour & Panfilov [TU] extend the saturation to 
artificial negative values, so that system Q does not 
degenerate in the one phase region and the saturation 
can still be used as a primary variable. 

J afire & Sboui [TT] use the solubility as an additional 
third primary variable. Additional nonlinear comple- 
mentarity constraints, which describe the transition 
from one phase to two phase region are used to close 
the system. 

We developed a new approach using the pressure of 
the nonwetting phase and the capillary pressure as pri- 
mary variables. In the absence of the nonwetting phase, 
p n is defined as the corresponding pressure to the sol- 
ubility x*. 

Our approach has the advantage, that we only have 
two primary variables in contrast to the complemen- 
tarity constraints method, where an additional variable 
is needed. With our constant set of variables we also 
avoid a switching of the primary variables, which is a 
non-differentiable process that can lead to numerical 
difficulties. 

This idea was first presented by Ippisch [13]. In the 
context of nuclear waste management for the special 
case that Henry's Law is used to couple solubility and 
pressure there exist similar approaches. Bourgeat et al. 
[3J use the water mass concentration and the wetting 
phase pressure, Angelini et al. [T5] use the two phase 
pressures as primary variables. 

In section [5] we will apply our approach not only to a 
recent benchmark study on nuclear waste management, 
but also to the very challenging field of CCS. In contrast 
to nuclear waste management and the work of Bourgeat 
et al. and Angelini et al. it is not possible to use Henry's 
Law for the solubility, because the approximation is 
not valid for CO2. We need a nonlinear function (see 



subsection 3.1 ) to describe the dependency between the 
nonwetting pressure and mole fraction. Moreover we 
have to handle the very high injection rate of the CO2. 



4.1 Pn/Pc formulation: Interpretation as algebraic 
transformation 

The entry pressure p on try is the critical capillary pres- 
sure that must be applied so that the nonwetting phase 
appears. We have to distinguish between 



to solve the problem at the phase transition (see Table 
[l] X£ denotes the mass fraction). 

A common method is primary variable switching 
used for example by Forsyth & Simpson fJ2] and 
Hclmig & Class [T3]. Here different sets of primary 



L Pc < Pcntry where S n — and only the wetting 
phase exists 

2. p c > p cri try where S n > and both wetting and 
nonwetting phase exist. 
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Case 1: 



Pc — Pcntry 



As mentioned in the beginning of the section the natu- 
ral set of variables for the one phase system ^ would 
be p w and x^. Consider the following transformation 
of variables 



Pw = Pn - Pc 



(10) 



where -0 is a continuous and invertible function. The 
solubility relation for x% ^ satisfies these demands 
(see Figure [I] and Spycher & Pruess [3]). The mapping 
between p n and x w is hence unique and p n and p c is a 
valid set of primary variables. 

The relation between the capillary pressure and the 
saturation p c (S w ) (see Equation ^) is a strictly de- 
creasing function for S w € [0, 1] and can therefore be 
inverted 

S w = viPc)- 

The dependent variables are then obtained through 





= rj(p c ) S n = 1 


- V(Pc) 


x h 


= V>0n) < = 1 


~ ^(Pn) 


X & 


= 7(Pn) = 1 


~ liPn) 



of the p n /p c formulation is, that the pressures, in con- 
trast to the saturations, are continuous across material 
heterogeneities if both phases exist. 

5 Numerical simulation 

In the following section we present the numerical results 
for special test cases. All simulations were performed 
in the DUNE simulation framework [TB], [T7] . 

A cell-centered finite volume method with two-point 
flux approximation on a structured grid was used for 
the domain discretization. The grid Eh — e%, . . . ,e n 
consists of elements e,; and the boundary of each ele- 
ment is dei = Ujgs(i) lij where 7^ denotes the bound- 
ary between elements and ej . The cell-centered finite 
volume method for Equation ^ for each component k 
then reads 

< / <f>d t {pmoi.w x w S w + p moi , n XnS n }de 

e t EE h l I 

H~ 7j iT ^ ^ / ' {Pmol,w *^w^w ~^ Pmol,n ^n^n} 

l|ei " jem i } 

■ V •{.,;;. • </')•/,,.. d 7 l =0, 



where 7 is the solubility curve given in Q. All other 
variables are computed as given in Section [3j 

This choice is not unique, another possible set would 
be p w /pc or p w /pn- Using p n as a primary variable has 
the advantage, that the highly nonlinear density and 
viscosity functions are directly dependent on a primary 
variable. We prefer p c over p\ as additional primary 
variable, because then the saturation only depends on 
the primary variable p c through the nonlinear capillary 
pressure-saturation relationship. 

Instead of the nonwetting phase pressure p n the mo- 
lar fraction x w — ip(p n ) could also be used as primary jg determined 
variable, which is very similar to the water mass con- 
centration used by Bourgeat et al. [5] . 



where n$j denotes the unit outer normal to 7^. 

A special upwinding scheme is used to calculate the 
phase fluxes at the interface between two elements to 
handle material discontinuities resulting in different 
capillary-pressure saturation curves and relative per- 
meability functions in both elements. 

The direction of the flux of phase a at the interface 
between two elements i and j can be obtained from the 
sum of the pressure gradient and the force of gravita- 
tion W a ^ij — {^Pa Pmass,a,ij ' <?) ' ^ij , where Pm&ss,a,ij 

is computed as the arithmetic average of cells and 
"j. Depending on the sign of w a ^j the upwind element 



upwind c 



Case 2: p c > p cntr y 

The common choice of primary variables in the two- 
phase region is one pressure and the saturation. 
With the Pn/Pc formulation we obtain the saturations 
through the retention curve S w = rj(p c ), the other vari- 
ables are computed accordingly. 

p w and x w are continuous at the interface between 
the one-phase and the two-phase region. Through the 
transformation (10), p n and p c are continuous at the 
interface too. With p n /p c we found a set of primary 
variables that can be consistently used in the presence 
or absence of the nonwetting phase. One advantage 




> 



The capillary pressure of the upwind element is used 
to calculate the relative permeability in each element. 
The obtained relative permeabilities are multiplied by 
the absolute permeabilities and the viscosities in each 
element. A harmonic average of the values is used to 
calculate the flux at the interface: 



K aii = K„ 



&ra.i (Pc,upwind a ) 



K ad = K j 



kra,j (Pc,upwind a ) 
K K 
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For homogeneous porous media this upwinding scheme 
corresponds to an upwinding of saturation. 

For the calculation of the convective component 
transport a full upwinding of the molar fractions and 
the molar densities based on the upwind direction is 
used with 



Value 



Value 



x 



a,upwind 



Pmol,a,ij — Pmol.a,upwind 



As time discretization scheme the implicit Euler 
Method was used. Newton's Method was applied to 
linearize the system. The Jacobian matrix is derived 
through numerical differentiation. The resulting lin- 
ear equation system is solved with a BiCGStab solver 
preconditioned by an algebraic mul-tigrid method (see 

BSD- 

We chose three different test cases, the first one is 
from a recent benchmark study concentrating on ap- 
pearance and disappearance of phases in the context of 
nuclear waste management. As there are no analytical 
solutions for two-phase two-component flow systems, 
we use the results of other groups as possibility to val- 
idate our results. We also conduct a grid convergence 
study to verify the experimental order of convergence 
of our implementation. 

With the second test case we apply our formulation 
to a CO2 sequestration scenario in 2D and perform a 
strong scalability test. The third test case extends the 
second test case to 3D and shows that our approach 
can handle the large number of unknowns. 

The simulations were performed in parallel with up 
to 16 processes. 

6 Test case 1: Gas injection in a fully water 
saturated domain (quasi-lD) 

The first test case is an example from the MoMas 
benchmark on multiphase flow in porous media |19| . 
[20] . We converted the descriptions to match the vari- 
ables used in this paper. 

In this case the considered nonwetting component is 
hydrogen and the wetting component is water. The sol- 
ubility of water in the nonwetting phase is neglected: 
x 1 ^ = 0. Hydrogen is injected into the left part of a 
rectangular domain (200 m x 20 m) with a flux of q" 1 
for 5 ■ 10 5 years. The domain is initially fully saturated 
by the water phase, consisting only of pure water with 
initial conditions p a — p a (see Table I2I). The boundary 
conditions are Neumann boundaries at the top and 
bottom (see Figure [H]). The Dirichlet boundary condi- 
tions for the outflow boundary are the same as the ini- 
tial conditions: p Q |p t — Pa- Gravitation is neglected, 
which leads to a quasi- ID problem. 

The relationship between p n and (where X* is 
the mass fraction in contrast to the molar fraction x 1 !) 
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Table 2: Parameters for test case 1 

nofiux 



noflux 



Figure 5: Test case 1: Domain setup 

is given through Henry's Law: 

H(T) 



X, 



Pi, 



where the partial pressure of hydrogen in the nonwet- 
ting phase is p^ = p n for this case because there is no 
water in the nonwetting phase for this example. The 
mass fraction X^ is then converted to the molar frac- 
tion: 

XlM a 



w X*M- + (1 - X*)M b ■ 

The nonwetting phase density is determined by the 
ideal gas law, wetting phase density is obtained through 
Henry's Law 



-'mass,!! 



y mass,w 



p n M h (RT)-\ 
The diffusion coefficient is given as 



10 3 + H(T)M n p b . 



pm,a 



4>s a 



M 3 



D K . 



A van Genuchten-Mualem model with the parame- 
ters n, a and i5 a)res as given in Table [2] is used for the 
soil water characteristic and relative permeabilities. All 
other parameters used in the simulation are also noted 
in Table [2 

A structured grid with 400 x 20 cells was used for 
the computations. Figure [6] and [7] show the nonwetting 
phase saturation and phase pressures at Ti n over time. 
S n is zero at the beginning, all injected hydrogen dis- 
solves into the wetting phase and no nonwetting phase 
is present. At t ~ 13000 years a nonwetting phase 
starts to appear at the injection point Ti n . 
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Figure 6: Test case 1: Nonwetting phase saturation at 
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Figure 7: Test case 1: Phase pressures at Th, 



For the computations we used a constant time step 
dt = 1000 years. We also verified the robustness of our 
model by using larger time steps (dt — 5000 years). 

Six different groups including our group participated 
in this benchmark example, the results of all groups 
are presented in [2U]. The results of our simulation 
corresponds well to the results of the other groups. 

In addition we performed a grid convergence study. 
For the coarsest level (level 1) we use 12 x 1 cells. 
For each level we double the amount of grid cells in 
a;-direction, so we have 12x2* cells for level i. The so- 
lution on level 12 was used as a reference solution. The 
resulting experimental order of convergence (EOC) can 
then be computed through 



EOC l+ i 



1 



log (2) 



log 



ei+i 



where is the L2-error between the solution on level 
i and the reference solution. At t = 2 ■ 10 5 years we 
get second order grid convergence for nonwetting phase 
pressure and capillary pressure (see Table [311 . 



level 


^elements 


EOC (p c ) 


EOC (p,,) 


1 


24 


2.01 


2.02 
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1.97 


1.98 
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1.98 


1.98 
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1.99 


•5 


384 


2.00 


1.99 


6 
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2.00 


2.00 


7 


1536 


2.01 


2.01 


8 


3072 


2.03 


2.02 



Table 3: Grid convergence study for test case 1 

The convergence study shows that our numerical so- 
lution converges with an optimal EOC of two, which is 
the maximum order that can be achieved with a cell- 
centered finite volume discretization. 

7 Test case 2: CO2 injection in a fully water 
saturated domain (2D) 

In the second test case CO2 is injected into the lower 
left part of a rectangular geometry (600 m x 100 m) with 
a flux of q™. The domain is located 800 m under the 
surface. As in test case 1, the top and bottom of the 
domain have noflux boundary conditions (see Figure 
[9l. For the Dirichlet boundary on the right side we 



noflux 



nonin 
Tin 



noflux 



Figure 9: Test case 2: Domain setup 

choose hydrostatic pressure for the water phase and 
zero pressure for the CO2 phase (which leads to x'F = 0) 

'• ' (900-z)p n 



P 



w| r„ 



= W 



y mass,w 



■g Pa, p h \ r =0Pa, 



where z is the z-coordinate in the domain and g the 
gravity in z-direction. Again the same values are taken 
as initial values. 
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Table 4: Parameters for test case 2 

Densities, viscosities and solubilities are chosen as 
suggested in Section [3j all other parameters are given 
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Figure 8: Test case 2: CO2 phase saturation and molar fraction of dissolved CO2 in water (contour lines 
for a& = 0.005, 0.011, 0.016). Color scale ranges from S n = (blue) to S n = max(5 , n ) (red). 



in Table [4] In this example we used the Brooks-Corey 
model for the soil water characteristic and relative per- 
meabilities. For the computations a structured grid 
with 240 x 40 cells was used. 

The results of test case 2 are shown in Figure [8] Each 
picture shows the CO2 phase saturation at a specific 
time point, the contour lines depict the molar fraction 
of CO2 in the water phase. The CO2 migrates upwards 
until it reaches the top of the domain with the noflux 
conditions and is then driven to the right by advective 
forces. Around the CO2 front the water phase contains 
dissolved C0 2 . 

An analytical solution for this test case does not ex- 
ist, but the simulation results are plausible and the 
C0 2 front behaves as expected. 

During the initial phase of CCS CO2 is injected and 
the CO2 front is propagating. Thus for the sake of 
accuracy we want to choose a time step size so that the 
CO2 front travels one grid cell layer per time step. 

For the computations we used a maximum time step 
of dt = 5000 s. The time step size is halved if the 
Newton solver did not converge, it is doubled until the 
maximum time step is reached in case of convergence. 
With this time step control we achieve an average time 
step size of dt = 3575 s and the CO2 front moves about 
one grid cell layer per time step, which fulfills the above 
condition. 

All simulations were done in parallel. To analyze the 
parallel performance of the simulations, we conduct a 
strong scalability test, where the global problem size 



^processes 


total time [s] 


efficiency 


1 


13975 


1 


2 


7763 


0.90 


4 


4151 


0.84 


8 


2658 


0.65 



Table 5: Strong scalability test for test case 2 



stays fixed and the number of processes is increased. 
The efficiency is defined as 



where T± is the time for the sequential method, p the 
number of processes and T p the time for the parallel 
method with p processes. Table [5] shows the results for 
test case 2 for a simulation time of 65 days. The total 
time needed for solving the problem scales very well 
with the number of processes. 

For a possible comparison with other implementa- 
tions we list some important performance indicators in 
Table [6] TS is the amount of time steps that were per- 
formed (successful and unsuccessful), for average and 
minimum time step sizes only the successful time steps 
were regarded. NI is the average number of Newton 
iterations per time step (successful and unsuccessful), 
where 10 NI are the maximum number of iterations 
that were allowed. Table [6] shows, that the average 
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(c) 14 days, (d) 18 days, 

max(Sn) = 0.82, max(i^) = 0.021 max(S n ) = 0.82, max(i^) = 0.021 

Figure 10: Test case 3: CO2 phase saturation and molar fraction of dissolved CO2 in water (contour lines 
for i» = 0.005, 0.010, 0.016). Color scale ranges from S n = (blue) to S n = max(S n ) (red). 



^processes 


TS 


av. dt 


min. dt 


av. NI 


1 


2249 


3579.7 


156.25 


3.9 


2 


2276 


3527.6 


156.25 


3.9 


4 


2205 


3593.1 


312.5 


3.9 


8 


2282 


3514.4 


312.5 


3.9 



Table 6: Average time step size and number of Newton 
iterations of the strong scalability test for test 
case 2 



time step size and number of Newton iterations stay 
almost constant for different number of processes. 



8 Test case 3: CO2 injection in a fully water 
saturated domain (3D) 

For test case 3, we use the same parameters and a very 
similar setup as in test case 2 in Section [7j The differ- 
ence is that we look at a 3D domain as shown in Figure 

HU 

The domain is a cube with dimensions 100 m x 
100 m x 100 m. For the computations a structured grid 
with 60 x 60 x 60 cells was used. 

The results of test case 3 are shown in Figure [T0| 
As in test case 2 each picture shows the C0 2 phase 
saturation and the solubility of CO2 in the water phase. 
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Figure 11: Test case 3: Domain setup 
9 Conclusion 

We suggest a new method to deal with the problem 
of disappearing nonwetting phase in two-phase two- 
component flow simulations. We use the nonwetting 
phase pressure and capillary pressure as primary vari- 
ables. This allows us to use the same variables for both 
the monophasic and diphasic case, no switching of pri- 
mary variables is needed to treat the nonwetting phase 
appearance problem. For the special case of CSS, we 
specify our choices for the necessary constitutive rela- 
tions. 

We confirm our new choice of primary variables with 
numerical simulations for different test cases in 2D and 
3D. We simulate the special case of CO2 injection in ge- 
ological formations and took part in the MoMas bench- 
mark on multiphase flow, where hydrogen flow in nu- 
clear waste repositories was examined. All simulations 
are performed in parallel and scale very well with the 
number of processes. In the benchmark case our output 
corresponds very well to the results of other groups. 

Next we want to extend our simulations to a non- 
isothermal model and use adaptive grid refinement. We 
want to use massive parallel computing in order to sim- 
ulate realistic CSS scenarios with very large domains 
and long time spans. 
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